Internal state of granular assemblies near random close packing 
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The structure of random sphere packings in mechanical equilibrium in prescribed stress states, as 
studied by molecular dynamics simulations, strongly depends on the assembling procedure. Friction- 
less packings in the limit of low pressure are devoid of dilatancy, and consequently share the same 
random close packing density, but exhibit fabric anisotropy related to stress anisotropy. Efficient 
compaction methods can be viewed as routes to circumvent the influence of friction. Simulations 
designed to resemble two such procedures, lubrication and vibration (or "tapping") show that the 
resulting granular structures differ, the less dense one having, remarkably, the larger coordination 
number. Density, coordination number and fabric can thus vary independently. Calculations of 
elastic moduli and comparisons with experimental results suggest that measurable elastic properties 
provide information on those important internal state variables. 

PACS numbers: 45.70.-n, 45.70.Cc, 71.55.Jv, 81.40.Jj 



The mechanical properties of solidlike granular pack- 
ings and their microscopic, grain-level origins are an ac- 
tive field of research in material science and condensed- 
matter physics P, 13 , with practical motivations in soil 
mechanics and material processing, as well as theoretical 
ones, as general approaches to the rheology of "jammed" 
systems [2| are attempted. For long, the only accessible 
information on the internal state of granular packings has 
been the density or the solid volume fraction (f>. Its im- 
portance was recognized, e.g. in soil mechanics 0,0, and 
for the elaboration of concrete mixes |^. The concept 
of random close packing (RCP) is deemed relevant 
in many contexts. For monodisperse spheres, the RCP 
volume fraction value ^rcp — 0.637 was consistently re- 
ported 0,0- Discrete numerical simulation proved a 
valuable tool to investigate the internal state of packings, 
as it is capable to reproduce mechanical behaviors 
and to identify relevant variables other than $, such as 
coordination number and fabric (or distribution of con- 
tact orientations) 0, . The influence of the sample 
preparation procedure on the mechanics of solid granu- 
lates is widely recognized as crucial, in numerical exper- 
iments as well as in real ones, but it is seldom discussed 
in detail. Experimentally, a variety of techniques such as 
controlled rain deposition under gravity , or layerwise 
tamping are used. Numerically, loose contactless con- 
figurations (^^CTanular gases") are homogeneously com- 
pressed J^l llJI or dropped in a container under their 
weight jlSlllq. In both cases, configurations are usually 
classified by their density, other state parameters being 
implicitly regarded as density-dependent. The present 
Letter reports on a molecular dynamics study that fur- 
ther investigates those initial states and their dependence 
on the assembling procedure in the case of monodisperse 
spherical particles. Its aim is twofold: first to revisit (af- 
ter others 0|) the prevailing notion of "random close 
packing"; then to study the internal states obtained by 
procedures designed to approach such RCP states, and 



to suggest possible means to measure important hidden 
information on their geometry. 

A foreword is necessary about the mechanical defini- 
tion [l^ of a configuration of maximum density. Let us 
consider a collection of a large number of identical rigid 
spheres, and enclose them in a rectangular parallelipi- 
pedic box. The degrees of freedom are the positions of 
sphere centers and the lengths L" {a — 1, 2, 3) of the 
container edges parallel to coordinate axes (changes of 
cell shape might also be considered). Using arbitrary 
reference lengths Lq strain parameters = ^"l"" 
defined. A local minimum of volume V, or maximum of 

is obtained on minimizing, under impenetrability con- 
straints, the potential energy W = ~PJ2a that corre- 
sponds to external pressure P > 0. Constraints entail the 
definition of repulsive, normal contact forces as Lagrange 
parameters. Local maxima of $ in configuration space 
(called "strictly jammed states" in ^?|) are thus equiva- 
lently characterized as stable equilibrium configurations 
of rigid, frictionless particles under an isotropic pressure. 
This duality between contact forces and rigid constraints, 
presented e.g., in jl8l |. leads to define any such configura- 
tion, devoid of crystal nucleus (wh ich might be checked 
with suitable order parameters [19|), as a random close- 
packed state. Not surprisingly, it is common practice, in 
order to obtain dense states easily, to set friction coeffi- 
cients to zero in molecular dynamics calculations 0, ^3 ■ 

Real materials are made of elastic, deformable grains, 
a feature most often taken into account in simula- 
tions 0, 0, 0|. In those presented below, we imple- 
ment the Hertz law [i^ for normal forces, and, in the 
presence of friction, like in refs. 0, 0, 0| , a simplified 
form of the Cattaneo-Mindlin-Deresiewicz relations ^3 
(as suggested in for tangential forces. Stable equi- 
libria then minimize a potential energy containing an ad- 
ditional elastic term. Solutions for perfectly rigid grains 
are recovered in the limit of large contact stiffness. 
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relative to the level of externally imposed forces. Specif- 
ically, the typical ratio h/a of elastic normal deflection 
of surfaces in contact h to sphere diameter a is of order 
(P/Ef/^ for Hertzian contacts (force F ~ Ea^/'^h^/'^). 
Our simulations, to allow for comparisons with experi- 
ments, correspond to glass beads with Young modulus 
E = 70kPa and Poisson coefficient v = 0.3. The method 
is similar to that of refs. ^ (21, 0| except that stresses, 
rather than strain rates, are controlled (as in ref. [j^V 
Unless otherwise specified, results are averaged over 5 
samples of 4000 beads enclosed in a periodic cubic cell. 

Compression of a frictionless granular gas to equilib- 
rium under P = lOfcPa yields $ = 0.637±0.002 ~ ^rcp, 
like in refs. 0, H^. The difference with the value <&* 
one would obtain in the P ^ limit might be estimated 
to fisrt order in the small elastic deffections he of the 
contacts, via the theorem of virtual work [l^ . on equat- 
ing the work of the external pressure, PV{^ — $*)/<!> to 
the work of normal contact forces /c, X^c/c^c, between 
equilibrium configurations at P = and P = lOfcPa. 
One finds $ - $* ~ 1.15 x lO^'*, less than statisti- 
cal uncertainties on $. Configurations in the P — s- 
limit are endowed with particular properties related to 
absence of force indeterminacy and stability, which to- 
gether entail the isostaticity of the force-carrying struc- 
ture , which consequently has a coordination number 
z* = 6. On the force network at lOfcPa each ball has 
an average of z* = 6.073 ± 0.004 contacts, in agreement 
with refs. 0, 0, H^- This value excludes "rattlers", 
spheres that carry no force, which represent a fraction 
/o = (1.3 ± 0.3)% of the total number, and differs from 
the global coordination number, z = z*(l — /o) (obtained 
on attributing zero contact to rattlers). In all cases con- 
sidered, only a few ( 0.3%) isolated particles belong to 
crystalline regions, as defined in ref. |l9j, p. 4. 

Numerical configurations made with friction are known 
to exhibit lower densities and coordination numbers 0, 
Q, l23j |. In practice, compaction strategies should 
therefore avoid the mobilization of friction in contacts. 
Our definition of RCP states is naturally associated with 
an isotropic pressure. What about other possible states 
of stress ? To investigate those, isotropic samples of 
frictionless were submitted to stepwise axisymmetric tri- 
axial compression , i.e., keeping the coordinate axes as 
principal stress directions, (73 was increased, at constant 
P = (cti +(J2 -l-cr3)/3, by steps of 0.02P, while maintain- 
ing CT2 = (Ji- At each stress step one waits for mechan- 
ical equilibrium (to accuracy lO^^'a^P on the net force 
on each particle). With perfectly rigid balls, one would 
have a local minimum of = — (Ja^a- The maxi- 
mum principal stress ratio supported by frictionless balls 
was observed to reach 1.2 to 1.24. For each equilibrium, 
z* , /o, strains Cq (taking the initial isotropic state as 
reference), and the fabric parameter % =< 3 cos^ 9 — 1> 
{9 is the angle of normals to contacts with the major 
principal stress direction z) are recorded. Figure Q dis- 
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FIG. 1: (Color online) (a) Axial strain £3 and (b) volume 
fraction $ vs. stress ratio a^/ai. Small symbols: earlier 
simulations, 1372 balls. Connected dots: 4000 ball samples. 



plays 63 (the axial strain) and $ as functions of a^/ai. 
Stress-strain relations are elusive in those frictionless sys- 
tems, as axial strain responds quite erratically to deviator 
stress increments (the strain/stress curve was interpreted 
as the trajectory of a Levy flight in 2D [31). $ varia- 
tions with deviator stress are extremely small, with a 
slight increasing trend, so that observed densities remain 
within the small interval of reported RCP values. This 
lack of dilatancy, contrary to the classical Reynods [l^ 
argument, might appear as a paradox, as density should 
first increase on destabilizing density-minimizing config- 
urations. Such density increases are however likely to 
vanish in the large system limit. This lack of volume 
change agrees with experimental observations, as stresses 
in alleged RCP states {e.g., under gravity) cannot be al- 
ways isotropic. Likewise, the numerical simulations that 
should produce anisotropic stress states, due to a final 
oedometric compression 14j, or to gravity [isl. Ilij. nev- 
ertheless yield $ ~ ^rcp in the frictionless case. On 
the other hand, one observes (fig. a gradual build- 
ing of a moderate contact network anisotropy, appar- 
ently determined by stress anisotropy for macroscopic 
samples (a likely consequence of the "fragility" of equi- 
librium states [2JI). The ratio of vertical (parallel to axis 
3, C33) to horizontal (Cn) computed oedometric elastic 
moduli reach about 1.2 at the largest deviator. Conse- 
quently, assuming full success of the compaction method 
in avoiding friction mobilization, observed RCP configu- 
rations belong to a continuum of different states. Those 
differ at least by one internal variable: fabric anisotropy, 
related to stress anisotropy, which should be accessible 
via measurements of elastic constants. 

Let us now turn to systems in which the effects of 
friction are not entirely suppressed, and investigate, fo- 
cussing on isotropic states, the possible effects of two 
different compacting strategies. One is lubrication, ei- 
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FIG. 2: (Color online). Fabric parameter x versus principal 
stress ratio. Symbols as on fig.0 x is larger for contacts car- 
rying larger than average forces. Note the apparent regression 
of fluctuations as sample size increases, unlike on fig. 



ther with a viscous, fluid oil (hydrodynamic lubrication), 
or with a very thin, greasy sohd layer coating the grain 
surfaces. The suppression of friction is then expected 
to be more efficient in the course of the assembling pro- 
cess than in the static packing under higher loads. To 
produce numerical configurations bearing some similar- 
ity to experimental, lubricated ones, we simulated slow, 
isotropic compressions, with a high level of viscous dissi- 
pation, using a small friction coefficient /ig = 0.02 up to 
equilibrium states under P — lOkPa. To study the effect 
of higher pressures (maintaining isotropy), the friction 
coefficient was raised to /i = 0.3, a typical value for a 
dry contact. Such numerical configurations are referred 
to as B samples in the sequel, A samples denoting those 
that were assembled without friction (note, however, that 
elastic moduli were subsequently evaluated in the pres- 
ence of tangential forces, assuming type A configurations, 
well annealed at each pressure level, mobilize tangential 
elasticity in response to small load increments). 

Another common procedure to circumvent friction and 
produce dense samples consists in vibrating, shaking or 
"tapping" l2^. In practice, this is most often how ex- 
perimental RCP states are reached with dry beads, for 
which a certain amount of intergranular friction is un- 
avoidable. This might be understood, as the replacement 
of permanent contacts by repeated collisions precludes 
the mobilization of friction, and therefore destabilizes ar- 
rangements that were equilibrated thanks to tangential 
contact forces. A third series of numerical configurations 
(series C) was produced in order to mimic the internal 
state of vibrated samples, as follows. First, A samples 
at P = lOfcPa were dilated, scaling all coordinates by 
a factor 1.005, so that all contacts opened. Then parti- 
cles were attributed random velocities and left to collide 
without energy dissipation (like the hard sphere fluid of 



ref. [l9j) at constant volume. This shaking stage lasted 
until each particle had undergone an average of 50 colli- 
sions. Samples were finally isotropically compressed with 
fi — 0.3 and strongly dissipative dynamics, to equilibrium 
at P = lOfcPa. 

A comparison of low-pressure equilibrated configura- 
tions A, B and C reveals that the "imperfectly lubri- 
cated" samples B do not quite reach the RCP solid frac- 
tion, as $ = 0.627 ± 0.0002 is observed at lOfcPa. How- 
ever, the coordination number of active grains, z* = 
5.79 ± 0.007, was only slightly reduced and the pro- 
portion of rattlers, (1.7 ± 0.2)%, rather similar. On 
the other hand, "vibrated" samples C are denser, with 
<& — 0.635 ± 0.001 very close to ^rcp, but their coordi- 
nation number is considerably smaller, z* — 4.56 ± 0.03, 
with /o = (13.3 ± 0.6)%. Remarkably, this shows that 
density and coordination number can vary independently 
for the same material, without fabric anisotropy. In fact, 
z* values in the C case are as low as in loose samples 
(series D, z* = 4.62 ± 0.01, $ = 0.606 ± 0.002 at lOkPa) 
prepared from a cold "granular gas" by direct isotropic 
compression with /i = 0.3 (see fig. O. Compressions to 
higher pressures (keeping fi = 0.3) and calculations of 
elastic moduli enabled comparisons to experimental re- 
sults. On figure|2| z* , /o, <& and velocities of longitudinal 
and transverse sound velocities, respectively denoted as 
Vp and V5, are plotted versus P in the lOkPa to lOOMPa 
range. Vp and Vs are deduced from computed bulk mod- 

uh (P) and shear moduh (G) as Vp = ^J{B + ^G)lp 

and Vs — \/ Gj p {p is the mass density of the packing) . 
Sound speeds are compared to experimental measure- 
ments on glass beads by Domenico jl^l, and by Jia and 
Mills li^. Numerical results with C (vibrated) samples 
are clearly in better agreement with experimental results 
on dense packings of dry beads, than A samples(i.e., ini- 
tially prepared in an ideal isotropic RCP state) , for which 
Vp and Vs are too large in the lOOfcPa range, and vary 
too slowly with P. (Results for A samples are in excellent 
agreement with ref. [T^). Sound speed dependence on P 
can be approximated with a power law, the exponent be- 
ing roughly 1/6 for A samples, as predicted by simple 
"effective medium" approaches whereas it is close 
to 0.25 for Vs and about 0.22 for Vp in C samples (due 
to faster variations of modulus G). Despite the remain- 
ing uncertainties and discrepancies in the comparison to 
experiments (the numerical procedure is a caricature of 
the experimental assembling process, and no information 
is available on sample anisotropics in '28'!, while nu- 
merical configurations are isotropic) we conclude that the 
numerical simulations of dense samples of dry particles 
should probably involve such annealing stages, which re- 
sult in much smaller coordination numbers. 

Interestingly, Jia et al. ^\ also prepared packings with 
a very small amount of a solid lubricant (trioleine) and 
observed smaller densities than for dry packings (typi- 
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density. Solid fraction, coordination number and fabric 
can vary independently. Comparisons with laboratory 
measurements of sound speeds under varying pressure on 
glass bead samples show good correspondence with real 
assembling procedures, and suggest that elastic proper- 
ties provide access to coordination and fabric. 
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FIG. 3: (Color online) (a) Solid fraction (b) coordination 
number z*j(c) speeds of sound Vp and Vs (with experimen- 
tal points |27l 123) for samples of types A, B, C, D, versus 
isotropic pressure P. Vp and Vs correlate to z* rather than 
to Error bars are smaller than symbols. 



cally $ ~ 0.625 instead of 0.64 at lOOFa), but larger 
sound speeds (fig. OJ- Our results for B samples ("im- 
perfect lubrication") are qualitatively similar, as we find 
sound speeds close to A (ideal RCP) values, larger than C 
("vibrated") ones. Note also the agreement for the slope 
of log (T^) versus log(P) (close to 1/6). The larger elastic 
moduli of initially lubricated samples might thus be at- 
tributed to their higher coordination number (in spite of 
their lower density). 

To conclude, we propose to define ideal random close 
packings of spherical balls as equilibrium states devoid of 
crystal nuclei that remain stable without friction, and to 
regard compaction procedures as recipes to suppress fric- 
tion mobilization. Simulations reveal that, even though 
they all share the same volume fraction and coordina- 
tion number (in the small stress limit), such RCP states 
differ by at least another state variable, a fabric ten- 
sor: they develop anisotropic contact networks to sup- 
port anisotropic stress states. Simple modelling of two 
imperfect compacting strategies, vibration and lubrica- 
tion, show that isotropic configurations with very dif- 
ferent coordination numbers coexist close to the RCP 
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